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ABSTRACT 


The electron impact ionization and dissociative ionization cross section 
data of CO, C0 2 , CH±, NH 3 , and S0 2 , measured in our laboratory, have 
been parametrized utilizing an empirical formula based on the Born approx- 
imation. For this purpose a x 2 minimization technique was employed which 
provided an excellent fit to the experimental data. 
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I. INTRODUCTION 


Cross sections for the production of positive ions by electron impact on 
molecules find application in a wide variety of plasmas. The collision process 
may be represented by the following relations: 


e~ + MN MN + + 2e~ (1) 

M * + N + 2e~ (2a) 

4M + JV+ + 2e~ (26) 

M + + N+ + 3<“ (3) 

■Z* Species of higher stages of ionization, (4) 


where a is the cross section, MN is a diatomic molecule with M and N 
as component atoms and the “+” sign indicates a positive ion. Although 
equations (1) through (4) have been written for a diatomic molecule, similar 
relations hold for polyatomic molecules. 

There are various definitions of the cross section, a. They are given 
below: 

a) Partial cross section, <r p : It represents the individual process given by 
equations (l) through (4). 

b) Total cross section, a(T) : The total cross section is defined according 
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to the method used for its determination: 


i) If the particle counting method is employed for obtaining the total 
number of ions (singly as well as multiply ionized) produced as a result of 
electron impact then it is known as total counting ionization cross section , 
cr c (T), and is given by: 


<Mr) = X> p + X>*, (5) 

P » 

where cr p is the partial ionization process described by eqs.(l) through (3) 
and <r‘ is the partial cross section for the i th stage of ionization. 

ii) If the cross section data are generated by measuring the total ion 
current then the total ionization cross section, cri(T), is as follows: 


c ! (T) = Y,°P + Y, Z < a v' 

P i 


( 6 ) 


where is the stage of ionization. 

Cross sections for ionization have been measured since the 1930’s. The 
various methods employed, in the past, for this purpose have been described 
in detail in several review articles 1,2,3 ’ 4 previously published in the literature. 
Theoretical calculations for these cross sections are difficult due to many 
channels in the continuum contributing to the ionization process. Recently, 
however, the status of the theory of ionization of atoms and ions has been 
reviewed by Younger 6 . A similar survey for molecules has also been made 
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by Rudge 6 . 


The first approach to the calculation of ionization cross sections was 
based on the classical theory and was presented by Thomson 7 in 1912. His 
expression for the ionization is as follows: 

ffp = 4 *(y-)V l (l ( 7 ) 

where £ is the number of electrons in the target with binding energy I, X « 
El * 1 is the reduced ionization energy and I H is the ionization energy of the 
hydrogen atom. 

After Thomson, several classical, semi-classical, and empirical formulas 
have been proposed for the calculation of ionization cross sections. They 
are by: Elewert ( 1952 ) 8 , Gryzinski ( 1959 )®, Post ( 1961 ) 10 , Drawin ( 1961 ) U , 
Burgess ( 1963 ) 12 , Stabler ( 1964 ) 13 , Seaton ( 1964 ) 14 , Gryzinski I ( 1965 , Sim- 
ple ionization) 16 , Gryzinski ( 1965 , Double ionization) 16 , Vriens ( 1966 ) 17 , Lotz 
( 1967 ) 18 , McFarland ( 1967 ) 1 ®, Jain and Khare ( 1976 ) 20 , Green and Sawada 
( 1972) 21 and Bell et al. ( 1982 ) 22 . 

Quantum mechanical calculations by Bethe 23 , however, showed that the 
simple asymptotic behavior of <r p (eq.7) is incorrect and it should vary as 
*2^ at high electron impact energies. Based on Bethe’s theory, Bell et al. 22 
proposed an empirical formula for atoms and ions to fit the experimental 
data at all energies of the colliding electron. It is as follows: 
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*«-hH?) + i>h)1 ,8 > 

where A and o< are fitting coefficients and all other quantities have been 

defined previously. The above formula takes care of the behavior of a p {E) at 
high electron impact energies. The coefficient A can be calculated by fitting 
to the Bethe relation at high energies: 

a p (2?)=~[Aln(£) + S]- (9) 

It can also be obtained from the following relation: 

oo 

A=-f^-dE (10) 

na J & 

I 

where a P h is the photon-ionization cross sections for CO, CO 2 , CH 4, NH 3 , and 
S0 2 . For the sake of convenience of the modelers to utilize our data, we have 
parametrized them using equations (8) and (9). Experimental apparatus and 
procedures for obtaining these data are described in appendix I. Section II 
describes the fitting procedure, and in section III, results are presented and 
discussed. 
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II. PROCEDURE FOR PARAMETRIZATION 


a) Methods 

Cross sections for electron impact ionization of CO, C0 2 , CH 4 , NH 3 , and 
S0 2 were previously measured in our laboratory. Equation (8) was then 
fitted to these cross sections by a x 2 minimization technique. The details 
are described below. For in-depth discussion of the methodology, refer to 
Bevington 24 . 

For this purpose, we define a parameter x 2 in the following way: 

X 2 = ~ M ~2 !/«(*♦•) - /«(*»)] 2 > ( n ) 

t * 

where M is the number of parameters, N the number of data points, ^ the 

weight for each data point, /„(*,•) the experimental data at and /*(*,-) is the 
theoretical value at x<, calculated from eq.(8). Since x^’s (experimental data) 
do not change in the fitting procedure, / t (x) is a function of o,’s alone, i.e., 

/t(x) = / t (ax,a 2 a M )- Consequently, x 2 is also a function of the parameters 

a<’s alone, 


X 2 = X 2 (fll» a 2)' • • ) a \{)- 


5 



Thus we can obtain a minimal value for x 2 s at least a locally minimal 
value, by manipulating a<’s on the a { -space. There are many ways to minimize 
X 2 , two of which are used in this report. They are the gradient search and the 
linearization of the fitting function. These methods are described in detail 
below: 

i) Gradient search 

In the gradient search, the negative gradient of x 2 

M 

-vxi = J2 6a < s <> ( 12 ) 

i 

is calculated at some point (of, of a° M ). The parameters a<’s are then incre- 

mented simultaneously by 5o<, that is in the direction of the negative gradient 
-Vx 2 . Si, in this case, is the unit orthogonal vector in the direction of o< in 
the oi-space, 5o< is the a< -component of the gradient -Vx 2 . 

Expanding x 2 > using Taylor’s series expansion, as a function of parame- 
ters o<’s, we have 



( 13 ) 


The optimum values for parameter increments, 6a/ s, are those for which x 2 
is a minimum, that is when 
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^! = M + y' 

dak do-k 


We have, from eq.(14), a set of M simultaneous linear equations in 8a/ s. 
To solve for 8a/ s, we first express these equations in matrix form. Letting 


d 2 xl 


[ dajda k 


8a : - 


= 0, 


k=l...M 


(14) 



V 



6a = 



a = 


dai dai 


daida^f 


da M dai 


°lxl 


da^t daj>4 


and substituting eq.(l5) into eq.(l4), we have 


(15) 


ft = 6 a • a. (16) 

Multiplying both sides of eq.(16) by a -1 we have a matrix of the parameter 
increments 6a : 

6a = P • a -1 . (17) 
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ii) Linearization of the fitting function. 

Instead of expanding x 2 as in the gradient search method, we expand 
f t {x) as a function of parameters a,’s, 


M r 


ft{x) = fo(x) + J2 


d JM 8aj 

da,- 1 


( 18 ) 


Together with eq.(ll), x 2 can then be expressed as, 


5735 E £[/•<«> 


1 All A r 

= JfZJf £ ^ I /e(l<) “ fo{xi) ~ £ 

i * I i=i 


d fo(x) 
L da 3 


Sdj 


( 19 ) 


Minimizing x 2 with respect to 6a k , we have 


^x 2 

d(6a k ) 


-2 

N-M 



M r 


fe{ x i) fo{ x i) yi 


dfoM 


5 ay 


3/ 0 (xQ 

da k 


= 0, 


k=l...M (20) 


As before, we have £ = 5a • a, in matrix form, or in a more explicit form, 


- pi) 
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and 


N 


a ik « Y 


' 1 dfojxj) djofa) ' 

tr? day da k 


j, k = 1..M. 


( 22 ) 


b) Advantages and disadvantages of the above two methods 

Theoretically, the gradient search will take the search to the minimal 
point of x 2 - However, from a computational point of view, it presents 
some difficulties. Computing the gradient of x 2 requires the computer to 
go through many calculations. Thus at a near-minimal point, where the in- 
crements 6a } - ’s are considerably small, the gradient search has to be utilized 
many times and thus renders the method inefficient. The linearization of 
the fitting function, on the other hand, gives good results near the minimum 
point. However, by the nature of the Taylor’s expansion, this method is 
unreliable when x 2 is too far away from the minimum point. 

c) Gradient-Expansion search method 

It is obvious then that these two methods are complementary. As dis- 
cussed in the last section, the gradient search is effective when x 2 is at a point 
far away from the minimum point, and the linearization method is effective 
nearby. The best features of these two methods can be incorporated into one 
by increasing the diagonal terms of matrix a by a factor of A, so that: 

j5=fa-tt, (23) 
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where 


<*/k(l + A), if j k, ( 24 } 

a 3 k, j±k. ( ’ 

If A is small, the above equation is similar to that obtained from the 
linearization method. If A is large, the diagonal terms dominate the matrix, 
thus we have 

Pi « A • 6&i ■ oca (25) 

which yields similar solutions to those obtained from the gradient search 
method. 

d) Algorithm 

The above search method can be summarized as follows: 

1) Set A = 0.0001; Compute x 2 (a), 

2) Compute 5 a and x 2 (a + 5a), 

3) If x 2 (a+ 5a) > x 2 (a), increase X by a factor of 10 and repeat step 2, 

4) If x 2 (a + 5a) < x 2 (a), decrease X by a factor of 10 and return to repeat 
step 2, substituting a + 5a for a, 



where a = (oi, 02 , ... , a^), and 5a = (5oi, 502 , • • • , Sqm) 

The above steps can be terminated when the difference of x 2 (a + 5a) and 
x 2 (a) is smaller than some predetermined constant; 0.0001 was used in this 
report. 
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The programs, using the methods discussed above, were written in Pascal 
language and were executed on the HP9836C and Apple He computers and 
implemented by Caltech’s CS10 graphics library. The program for HP9836C 
is listed in appendix II. Flow diagrams from these programs are given in Fig- 
ures 1 and 2. 
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END 

Figure 2. Subprocedure CURFIT 
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m. RESULTS AND DISCUSSION 


Utilizing eq.(8) and the fitting procedure described in section II, the 
values of coefficients A and a,’s were determined. They are presented in 
Table I. Fitted data are also shown in graphical forms in figures 2 through 
26. Except for S0 2 , the experimental measurements were up to 510 eV 
electron impact energy. The data uncertainty was about 15%. The error 
bars are shown in all figures. In the cases of CO, C0 2 , CH i , and NH 3 , eq.(8) 
has provided an excellent fit to the experimental data. However, the data 
for S0 2 are only up to 200 eV and the fits are not very satisfactory. 
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TABLE I 


The parameters A and o<’s of (8) for the fitted cross sections. These 
parameters are in the units of 10 ~ 14 eV 2 cm 2 . 




A 

a x 

a2 

a 3 

a 4 

a 5 

a 6 


W2SM 

65.459 

-65.00 

-22.427 

-26.250 





W3SM 




6.189 

-6.089 



co 2 

c + 

1.433 

-1.262 

-2.514 

6.436 

-3.055 





9.355 

-9.252 

-5.309 

■KMA | 





EZ3SQI 

77.266 1 



27.582 





fgga 

49.342 

mmm 

-39.366 

45.765 

-51.574 



CO 

c+ 

6.781 


5.229 

11M 





o+ 

2.384 

-2.368 

0.984 

-10.776 

29.389 

-18.058 



Total 

52.961 

-51.003 

-35.790 

24.074 

-31.205 




waum 

25.660 

-24.545 

-9.057 

-1.971 





wmm 

30.133 

-27.286 

-43.680 

169.936 

-276.50 

140.548 


CH4 

W3ZM 

3.428 

-2.609 

-7.047 






WOTiM 

0.819 

0.0684 

-12.478 

88.954 

-219.99 

245.070 

-100.40 


C+ 

-0.0435 

0.330 

-1.828 


-74.772 




Total 

60.704 

-54.546 

-82.229 

287.768 

-447.52 

225.250 



wwm 

24.677 

-23.387 


-35.344 

65.208 

-48.343 



mum 

57.455 

-58.280 

-1.297 

-31.188 




EZ39 

NH+ 

4.780 




31.153 

-19.699 



N+ 

2.359 

-1.796 

-10.851 

68.642 

-186.8 

233.891 

-104.63 


Total 

60.008 

-59.746 

-15.417 

-88.234 

113.977 

■a riKatCT 




66.607 

-63.735 

-75.659 

331.36 





S0+ 

-9.180 

11.802 


403.755 




so 2 

5+ 

29.372 

-24.266 

-57.076 

196.126 


173.684 



o + 

2.010 


-2.537 

9.661 





Total 

80.341 

-77.366 

5.886 

-311.92 

1414.13 

-2309.7 

1252.17 
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Figure 3. Dissociative ionization and attachment spectrometer (dimensions not to scale) 
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Figure 4. Total electron impact ionization cross section for CO 
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APPENDIX I 

EXPERIMENTAL APPARATUS AND METHOD 


A schematic diagram of the apparatus is shown in Fig. 3. It utilizes a 
crossed electron beam — molecular/atomic beam collision geometry. In the 
cause of gases, the beam of atoms or molecules is produced by flowing the 
gas through a capillary array. Alternatively, in the case of species which 
are solids or liquids at room temperature, an electron bombarded oven or 
a resistance heated oven is utilized to produce the target beam. The beam 
of electrons is generated by heating a pure tungsten filament. The electrons 
are first extracted from the filament and are then accelerated or decelerated 
by three cylindrical lenses. This beam is collimated by the help of an axial 
B field which is produced by a solenoid within which the electron gun and 
a Faraday cup are housed. The solenoid produces the B field of the order 
200 G. 

The present beam of electrons is energy unselected. An energy profile of 
this beam was obtained by utilizing the retarding potential on the Faraday 
cup. It is found that the full width at half maximum (FWHM) is approxi- 
mately 300 meV. The energy of the electrons is varied by changing the bias 
on the filament with respect to the last electrode of the electron gun. This 
electrode is kept at ground (earth) potential. It was found that the beam 
current, as measured by the Faraday cup, remained constant as the energy 
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of the beam was changed from 0.5 to about 10 eV which is the range of the 
present interest. Although the energy of the electron beam can be obtained 
by recording the filament bias voltage, the contact potentials at various sur- 
faces tend to change it from its actual value. In the present work, the energy 
of the beam was calibrated by utilizing the accurately known values of ion- 
ization potentials of the rare gases. 

The positive ions produced by collision of electrons with the target mole- 
cules are extracted out of the B field by two parallel molybdenum wire meshes 
between which a voltage is applied. This voltage produces a homogeneous 
electric field with a gradient of 3 to 10 V/cm at the target. The direction 
of the field is normal to both electron beam and molecular beam. One grid 
is biased negative with respect to the ground and the other positive. The 
molecular beam is kept at ground potential. It was found that this arrange- 
ment did not disturb the electron beam. The efficiency of extraction of ions 
was measured by changing the extracting electric field from 0 up to 10 V/ cm. 
It was found that by increasing the electric field strength, the detected ion 
intensity increased rapidly in the beginning. However, at about 3 eV and 
above the ion intensity became almost constant as a function of the field 
strength. This indicated that the measured ion current did not depend on 
the initial energy and angular distribution of the ions. All our measurements 
were performed in this region of extracting voltage. The extracted ions are 
accelerated from 0 to about 200 V/cm and focused at the entrance aper- 
ture of a quadrupole mass spectrometer by an ion lens (Fig.l). This mass 
spectrometer has a resolution of approximately 1 amu. The mass analyzed 
ions are accelerated by a 3.2 kV potential and are detected by a spiraltron 
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multiplier. Each ion is counted as an event by a multichannel scaler. 

A vacuum of about 10~ 8 Torr was obtained when the gas forming the 
molecular beam was not flowed into the vacuum chamber. However, the 
pressure rose to about 10“ 7 Torr when the molecular beam was on. 


In order to obtain the absolute values of the cross sections, the relative 
flow 25 technique developed in our laboratory for collision cross section mea- 
surement was utilized. The method employs a measurement of the ratio 
of the intensity of the positive ions of the unknown species (for example, 
O + /S0 2 , S + /S0 2 , S0 + /S0 2 , SO% /S0 2 ) to that of a known species (for exam- 
ple, He + /He, Ne + /Ne, Ar + /Ar, or Kr + /Kr). At the same time, the flow rates 
and pressure behind the capillary array 25 are measured. First, the gas [ AB 
of Eq.(23) below] whose ionization cross section has to be measured is flowed 
through the capillary array and a beam is formed. The positive ion intensity 
I{B + ) is then measured. Subsequently, the gas AB is turned off and Ar is 
flowed through the capillary array. The positive ion current I{Ar + ) is again 
recorded. Providing that the measurement is performed under the conditions 
of molecular flow through the capillary array, the following relation is used 
to obtain the cross section: 



m(Ar) 1 1/2 / N(Ar) \ 

m(AB) J \N{AB)J ’ 


( 26 ) 


where m(AB) and m(.Ar) are molecular weight of respective gases, N(Ar) and 
N(AB) are the flow rates of the two gases through the capillary array, and K 
is a calibration constant which determines the transmission efficiency of the 


48 



ion optics, quadrupole mass spectrometer, and charged particle detector for 
B + and Ar + . 

The calibration constant K for the various masses was experimentally ob- 
tained. We chose gases whose ionization cross sections are well known. These 
are {H+/H 2 ) 26 , (He+/He) 27 , ( 0 + , 0 2 ) 26 , ( Ne+/Ne ) 27 , {Ar+/Ar) 27 , and Kr+/Kr) 27 . 
Since all the quantities in Eq.(23) are either known or can be obtained ex- 
perimentally for the two gases out of the ones mentioned above except K , 
the values of K for various mass numbers ranging from H to Kr can be calcu- 
lated. We followed this method for calibrating our instrument. The relative 
efficiency if as a function of the mass number is a bell shaped curve. It is 
increasing with the mass number up to about mass number 45, then it is de- 
creasing at higher mass numbers. Our results are in agreement with Ehlert’s 
measurements 28 . More detail about the relative efficiency measurement can 
be found in Orient and Srivastava’s paper 29 . 

The contribution of the background scattering (both direct beam con- 
tribution and scattering by the background gas) to the scattering from the 
target gas beam is measured by providing an alternate leak to the vacuum 
chambers. The flow to the chamber is switched from the capillary array to 
the alternate gas inlet and the proper background pressure for the desired 
gas is established. The mass selected beam intensity is then measured as a 
function of the electron beam energy. It is found that the maximum value 
of the background scattering is about 5%. 
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Program Listing for HP9836C Computer 
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program surf 1 ( input , output) ; 

import mylib; { graphic library C.I.T. } 

const 


top 

« 350; 

pi * 3.14159268; 

maxdata= 200; 

bottom 

= 20; 

maxpar = 30; 

termx = 9999 

left 

= 20; 

long * 480 ; 

maxfunc* 100; 

right 

= 500; 

height * 330; 

maxstack* 20; 

type 

chars 

- *#*. 

.'z*; 



charset = set of chars; 

elem.type = (cons , coef .oper.vars .ends) ; { type ol function elements } 

oper_type = (oadd , osub , omul , odiv , opwr , oexp , olog , oc os , osin , oatn , 

otan.ocsh.osnh.otnh) ; { defined basic operators > 

func_rec = record { function elements } 

case datum : elem_type of 

vars : (valr : real) ; 

coef .cons : (vali : integer) ; 

oper : (valo : oper_type) ; 

end; 


func_arr * array [0. .maxfunc] of func_rec ; { function array > 

stck_arr * array [1 . .maxstack] of real; 
vhtype = (vertical, horizontal) ; 

mode_type * (instr.none ,stat) ; { uncertainty type > 

theograf = array [left . .right] of real; { array of graph > 

coeff_rec * record { record of a coefficient > 

a.siga.dela : real; 
end; 

coeff_type= array [l..maxpar] of coeff_rec; { array of coefficients } 

cons_type * array [i..maxpar] of real; { array of function constants > 

data_rec * record { record of data > 

x.y.sigy : real; 
end; 

expdata * array [l..maxdata] of data_rec ; { array of data > 


filedata * file of data_rec ; 

f ilefuncrec»record { record of a saved function > 

sort , sub : integer ; 
end; 

filefunc * file of filefuncrec; 
mattype = array [1. .maxpar.l. .maxpar] of real; 
oneddata = array [l..maxdata] of real; 

var 

dmaxy, dminy, maxx, maxy, minx, miny : real; 
lamda,x2,xl ,corrx, corry, cutoff : real; 
datg , dat ; theograf ; 
indat, indatg : expdata; 
filea : filedata; 

cf terms , connum ,npts ,nf ree .color : integer ; 


mode : mode_type ; 
funca : func_arr; 
coeff ; coeff_type; 
name : string [12] ; 
fileb : filefunc; 
cona : cons_type; 
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mch : char; 


exitprog : boolean; 


{ > 

function ff(fop : oper_type ;xi ,xi_l :real):real; { elementary operations } 
begin { between two elements > 
case fop of { select operation > 


osub : ff : = xi_l - xi; 

oadd : ff : = xi_l + xi; 

omul : ff :■ xi_l * xi; 

odiv : ff :■ xi_l / xi; 

opwr : ff := exp(xi*ln(xi_l)) ; 

oexp : if abs(xi) < 700 then ff := exp(xi) else 

if xi > 0 then ff :« le300 else ff := le-300; 
olog : if xi > le-300 then ff :* ln(xi) else ff := -le300; 

ocos : ff := cos(xi); 

osin : ff sin(xi) ; 

oatn : ff :* arctan(xi) ; 

otan : ff :* sin(xi)/cos(xi) ; 

osnh : if abs(xi) < 700 then ff := 0.5*(exp(xi) - exp(-xi)) else 

if xi > 0 then ff le300 else ff := -le300; 

ocsh : if abs(xi) < 700 then ff := 0.5*(exp(xi) + exp(-xi)) else 

ff :* le300 ; 

otnh : if abs(xi) > 700 then ff := (exp(xi) - exp(-xi) )/(exp(xi) + exp(-xi)) 
else if xi > 0 then ff :■ 1 else ff : = -1; 

end; 

end; 


{ > 

function f (ffunca:func_arr;fx,f cutoff : real;fcoeff: coeff_type; 

fcona : cons_type): real; { evaluate function ffunca at > 

var i,st : integer; fstck : stck_arr; { fx > 

begin 

st :* 0; { stack count } 

i := 1; 

if fx O f cutoff then f := 0 else { f * 0 below cut off point > 

begin 

while ffunca[i] .datum <> ends do { unstack k evaluate function } 

begin 


if ffunca [i] .datum <> oper then st :« st + l;{stack unless operator} 
with ffunca [i] do 
case datum of 


vars 

f stck[st] 

• = 

fx; { 

load values onto stack 

> 

cons 

IstckCst] 

• s 

fcona[vali] ; 



coef 

lstck[st] 

; = 

f coeff [vali] 

.a; 


oper 

case valo 

of 


evaluate w/ given operator 

> 


oadd, osub, odiv, omul , opwr : begin 

f stck[st-l] :* ff (valo,f stck[st] ,f stck[st-l] ) ; 
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st :* st -1; 
end; 

oexp.olog.ocos , osin,oatn,otan,ocsh,osnh,otnh : 
fstck[st] :■ ff (valo ,f stck[st] ,f stck[st] ) ; 
end; 

end; 


i :■ i + 1; { increment stack count > 

end; 

f :* fstck[l] ; { f is first element of stack } 

end; { after evaluation > 

end; 

{ - > 

procedure maxmindata(mindat : expdata) ; { find maximum and minimum } 

var i : integer; { of both x and y of data > 

begin 

maxx :■ mindat[l].x; { initial max. min. } 

maxy :* mindat[l].y + mindat [1] . sigy; 

minx :* mindat [l].x; 

miny ■:* mindat [1] .y - mindat [1] .sigy; 

for i :* 1 to npts do { search for max. min. > 


with mindat [i] do 
begin 

if x > maxx then maxx :■ x; 

if x < minx then minx x; 

if y + sigy > maxy then maxy :■ y + sigy; 

if y - sigy < miny then miny :* y - sigy; 



end; 


dmaxy 

:■ maxy; 

{ save dmaxy, dminy t or graph > 

dminy 

:* miny; 


corrx 

: *= long/ (maxx - minx); 

{correction factors for screen} 

corry 

: * height/ (maxy - miny); 


end; 




{ - > 


{ - — -> 

procedure maxmin; { calculate max and min for > 
var i : integer; { both the graph > 
begin { and input data } 


maxy : * dmaxy ; 
miny : * dminy ; 
for i :■ left to right do 
begin 

if dat[i] > maxy then maxy :* dat[i] ; 
if dat[i] < miny then miny :■ dat[i] ; 
end; 

corry :* height/(maxy - miny); { recalulate y correction > 

end; 
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function theodata(tf unca : func_arr;tcoeff : coeff_type; 

tcona : cons_type) : theograf ; { caluculate values for graph > 

var i : integer; tcorx : real; 
begin 

tcorx :■ long/Cmaxx - minx); 

for i :■ left to right do { for every pixel points > 

theodata[i] := f(tf unca, minx + (i-left) /tcorx, cutof f .tcoeff .tcona) ; 
end; 


function theonorm(tdat : theograf): theograf; { normalize values for graph > 

var i : integer; 

begin 

for i := left to right do theonorm[i] := (tdat[i] - miny)*corry + bottom; 
end; 


function expnorm(eindat : expdata) : expdata; 

var i : integer; { normalize data for graphics > 

begin 

for i := 1 to npts do 
begin 

expnorm[i] .x := (eindat[i].x - minx)*corrx + left; 

expnorm[i].y := (eindat[i].y - miny) *corry + bottom; 

expnorm [i] . sigy := eindat [i] . sigy*corry ; 
end; 
end; 


procedure drawtheo(ddatg : theograf); 

var i : integer; 

begin 

m_move(left,round(ddatg[left])) ; 
for i := left to right do 
m_draw(i , round (ddatg[i] )) ; 
end; 


procedure drawexp(dindatg : expdata;dnpts : integer) ;{ draw data points } 

var i,dx,dy,dsigy : integer; 

begin 

for i := 1 to dnpts do 


{ draw the graph function > 
{ first point > 
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begin 

dx :* round (dindatg [i] .x) ; { rounding normalized values } 


dy :■ round (dindatg [i] .y) ; 
dsigy :* round (dindatg ti] . sigy) ; 

m_drawrect(dx-l , dy-1, dx+i, dy+i) ; { draw small box for point } 

m_move(dx-l, dy+dsigy) ; { draw end point for error bar} 

m_draw(dx+l, dy+dsigy); 

m_move(dx, dy+dsigy); ■( error bar } 

m_draw(dx, dy-dsigy) ; 

m_move(dx-l, dy-dsigy); •( another end point } 

m_draw(dx+l, dy-dsigy): 
end; 
end; 

i > 

procedure readdata(var rindat : expdata;var rnpts : integer ;strt: integer) ; 

var i : integer; { user input data } 

begin 

write ( ’Condition of sigma-y: 1-Statistical 2-No sigma-y 3 -Instrumental ’); 

readln(i) ; { kind of uncertainty } 

rnpts :* strt; { strt * 0 for new data } 


repeat 

rnpts :■ mpts+ 1; 
case i of 

1,2 : begin 

rindat [rnpts] , sigy :■ 0; 

if i * 1 then mode : * stat else mode : * none ; 
write(’ (’ .rnpts :1, ’) x, y? (x * 9999 to end) ’); 
readln (rindat [rnpts] . x , rindat [rnpts] . y) ; 
end; 

3 : begin 

mode : * instr ; 

write (* (’ .rnpts :1. *) x, y, sigma-y? (x * 9999 to end) ’); 
readln (rindat [rnpts] . x , rindat [rnpts] . y , rindat [rnpts] . sigy) ; 
end; 

otherwise ; 
end; 

until (rindat [rnpts] .x * termx) or (rnpts ■ maxdata) or (i<l) or (i>3) ; 


if rindat [rnpts] .x * termx then rnpts :* rnpts-i;{ do not store 9999 } 

end; 

i > 

procedure save! ile(sindat : expdata;snpts : integer); 

var i : integer ; { save input data } 

begin 

write ('name of file *); 
readln (name) ; 
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rewrite(f ilea, name) ; 
for i : s 1 to snpts do 
begin 

filea* :* sindat[i]; 
put(f ilea) ; 
end; 

close (lilea, ’save*) ; 
end; 

{ - - -} 

procedure readfile(var rindat : expdata;var rnpts : integer ;strt : integer); 
var i : integer; { read from file input-data > 

begin 

write ('read from '); 
readln(name) ; 
reset (f ilea, name) ; 
rnpts :* strt; 
while not eof(filea) do 
begin 

rnpts := rnpts+1; 
rindat [rnpts] := filea"; 
get (f ilea) ; 
end; 

close(filea) ; 

write(’Is sigma-y 1-Statistical 2-None 3- Instrumental *); 
readln(i) ; 
case i of 

1 : mode :* stat; 2 : mode none; 3 : mode instr; otherwise; end; 
end; 


i — - - — > 

function weight (windat : expdata;i: integer): real;{ calculate weights for > 

begin { each data points > 


with windat [i] do 
case mode of 

stat : if y <> 0 then weight := l/abs(y) 
else weight := 1; 
none : weight := 1; 

instr: if sigy <> 0 then weight :* l/sqr(sigy) else weight := leiO; 
end; 

end; 


{- — - - 

function chisqr (cindat : expdata; cfunc : func_arr; ccoef f : coef f _type) : real; 
var cx2 : real; •( calculate chi square } 

i : integer; 
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begin 
cx2 :* 0; 

ii nfree <* 0 then chisqr :* 0 else 
begin 

lor i:= 1 to npts do cx2 :* cx2 + weight (cindat ,i)*sqr(abs(c indat[i] .y - 


f (cfunc .cindat [i] .x, cut off , ccoeff ,cona))) ; 
chisqr :■ cx2/nfree; { reduced chi squared > 

end; 
end; 

{- - - - > 

function deriv(dindat : expdata;i : integer) : oneddata ; 

var j : integer; { calculate derivatives } 


coj ,yf it .durader : real; { with respect to all coeff. > 

begin 

for j := 1 to cf terms do 
begin 

coj := coeff [j] .a; 

coeff [j]. a : * coj + coeff [j] .dela; 

yf it :« f (funca.dindat [i] .x, cutoff .coeff ,cona) ; 

coeff [j]. a :* coj - coeff [j ]. dela; 

dumder : * (yf it-f (funca.dindat [i] .x, cutoff .coeff ,cona) )/(2*coeff [j] .dela) ; 
if dumder <> 0 then deriv[j] :* dumder else deriv[j] : * le-10; 
coeff [j] . a : * coj ; 
end; 
end; 

{ - — > 

procedure swap(var sl,s2 : real); { swap two numbers si and s2 } 

var dummy ; real; 

begin 

dummy : = si; 
si :* s2; 
s2 : = dummy ; 
end; 

{ > 

procedure matinv(var marray ; mattype;var merr : boolean); 

var i.j.k.l : integer; { invert matrice, calc. det. } 

arrmax. dummy : real; ik.jk : array [l..maxpar] of integer; 
begin 

merr :* false; 
for k :» 1 to cf terms do 
begin 

arrmax :■ 0; 

for i := k to cf terms do { find largest element > 
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{ swap column and row 
{ matrices to insure 
{ computational accuracy 


for j := k to cfterms do 

if abs (arrmax) <= absCmarray [i , j] ) then 
begin 

arrmax :* marray[i,j]; 
ik[k] := i; 
jk[k] := j; 
end; 

if arrmax = 0 then merr := true else 
begin 

i := ik[k] ; 
if i>k then 

for j :* 1 to cf terms do swap(marr ay [k, j] .marray [i , j] ) ; 
j := jk[k] ; 
if j >k then 

for i := 1 to cfterms do swap (marray [i ,k] .marray [i , j] ) ; 

end; 

if not merr then { invert matrix 

begin 

for i ;= 1 to cfterms do 

if i<>k then marray [i,k] := -marray [i ,k] /arrmax; 
for i := 1 to cfterms do 
for j := 1 to cfterms do 

if (iok) and (jOk) then 

marray[i,j] marray[i,j] + marray Ci. k] *marray [k, j] ; 
for j ;= 1 to cfterms do 

if j<>k then marray [k.j] := marray [k, j] /arrmax; 
marray [k,k] ; = 1/arrmax; 
end; 

end; 

if not merr then { return matrice to its 

for 1 := 1 to cfterms do { prearranged state 

begin 

k := cfterms -1+1; 
j : * ik [k] ; 
if j>k then 

for i := 1 to cfterms do swap (marray [i.k] .marray [i.j] ) ; 
i :* j k [k] ; 
if i>k then 

for j := 1 to cfterms do swap (marray [k, j] .marray [i , j] ) ; 

end; 

end; 


> 

> 

> 


> 

> 


-} 


procedure curf it (cindat : expdata;var ccoeff : coeff_type; 

var cx2 : real;cwt : oneddata) ; { calculate grad, fc new coeff .} 
var xsql : real; carray, alpha : mattype; 

cderiv : oneddata; i.j.k.l : integer; 

b : coeff_type; beta : array [1 . .maxpar] of real; 
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{ initialize matrice 


> 


cerr : boolean; 
begin 

lor j :« 1 to cf terms do 
begin 

beta[j] : m 0; 

lor k :« 1 to j do alpha[j,k]:» 0; 
end; 

lor i :* 1 to npts do { calculate beta and alpha > 

begin 

cderiv :■ deriv(cindat , i) ; 
lor j :* i to clterms do 
begin 

beta[j] :* beta[j] + cwt[i]* 

(cindat[i] . y-l(lunca,cindat[i] .x.cutoll , ccoell , cona) )*cderiv[j ] ; 
lor k :■ 1 to j do 

alpha [j ,k] :■ alpha [j ,k] + cwt[i]*cderiv[j]*cderiv[k] ; 

end; 

end; 

lor j :« 1 to clterms do 

lor k :»1 to j do alpha[k,j] :■ alpha [j ,k] ; 
xsql :« cx2; 

repeat { begin calculating direction } 

lor j :* 1 to clterms do { ol minimal chi square > 

begin 

lor k :* 1 to clterms do { normalize alpha lor accuracy} 

carray[j,k] :■ alpha[j ,k]/sqrt (alpha [j , j]*alpha[k,k]) ; 
carray[j,j] :* 1 + lamda; 
end; 

matinv (c array , cerr) ; { invert mtrix } 

il cerr then writeln( 'Non-unique system ol equation lound ’); 
b :* ccoell; 

lor j :■ 1 to clterms do { new coellicients? } 

lor k :* 1 to clterms do 
b [ j ] .a :■ b [ j ] .a + 

beta[k]*carray[j ,k]/sqrt (alpha [j , j]*alpha[k,k]) ; 
cx2 :« chisqr (cindat , funca.b) ; 


il cx2 > xsql then lamda :» lamda* 10 ; { gradient search? } 

until cx2 <* xsql; { keep going to minimum pt } 

ccoell :* b; { new coellicients } 

lor j :* 1 to clterms do ccoell [j] . siga :■ sqrtCcarray [j , j ] /alpha [j , j] ) ; 
lamda :« lamda/10; { linearization? } 

end; 

{ > 

procedure drawscale(dmin,dmax:real;dcon,dsc :integer;dvh:vhtype) ; 
var dlac : real; { draw scale on graph } 

pwcor , pwmax, pwmin, pwdil, gmin, gmax, i, k : integer; 

begin 
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if dmax > dmin then { compare magnitudes of max. 

begin { min. and calculate steps 

if abs(dmax) < 1 then pwcor := -1 else pwcor :* 0; { for scales 

if dmax <> 0 then pwmax : = trunc (ln(abs(dmax) ) /ln(10) ) + pwcor 
else pwmax :* -400; 

if dmin <> 0 then pwmin :* trunc (ln(abs(dmin) ) /ln(10) ) + pwcor 
else pwmin := -400; 

pwdif := trunc(ln(dmax-dmin)/ln(10)) + pwcor; 
if pwmax > pwmin then dfac :» exp((-l + pwmax) *ln( 10)) 
else if pwmin > pwmax then dfac :* exp((-l + pwmin) *ln(10)) 
else dfac :* exp((-l + pwdif ) *ln(10) ) ; 
gmax :» trunc (dmax/dfac ) ; 
gmin :* round (dmin/dfac + 0.5); 
for i :* gmin to gmax do 
begin 

case dvh of 

vertical : begin { left and right verticals 

k :■ round(corry*(i*dfac - dmin)); 
m_move(dcon,k + bottom); 
if (i*0) and (dcon ■ left) then 
m_draw (right ,k + bottom) else 
if i mod 10 * 0 then m_draw(dcon+3*dsc ,k+bottom) 
else m_draw(dcon+dsc ,k + bottom); 

end; 

horizontal: begin { top and bottom horizontals 

k :* round(corrx*(i*df ac - dmin)); 
m_move(k + left, dcon); 
if (i*0) and (dcon = bottom) then 
m_draw(k + left, top) else 
if i mod 10 ■ 0 then m_draw(k + left , dcon+3*dsc) 
else m_draw(k + left ,dcon+dsc) ; 

end; 


end; 

end; 

end; 

end; 


> 

> 

> 


> 


> 




> 


procedure drawframe; { draw rectangle w/ scale } 

begin 

m_color(m_yellow) ; 

m_drawrect(left .bottom, right , top) ; { draw box > 

writeln( ’maxy * *,maxy,’ miny = ’,miny); { draw scales > 

drawscale (miny , maxy , left , 5 , vertical) ; 

drawscale(miny, maxy .right , -5, vertical) ; 

drawscale (minx , maxx , top , -5 .horizontal) ; 

writeln( ’maxx * ’,maxx,’ minx * ’.minx); 

drawscale (minx .maxx , bottom , 5 .horizontal) ; 
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m_color(m_red) ; 

m_move (left , 10) ; 

m_di splay text ( ’Data file: ’); 

m_displaytext(name) ; 

end; 

{ > 


procedure draweverything(dindat : expdata; var ddat : theograf) ; 

var dindatg : expdata; ddatg : theograf; { draw both data and function } 

begin 

maxmindata(dindat) ; { calculating... } 

ddat :« theodata(funca,coeff , cona) ; 

maxmin; 

dindatg :» expnorn(dindat) ; 
ddatg : * theonorm(ddat) ; 

m.clear; < and drawing > 

drawf rame ; 

m_color (color) ; 

dr awtheo (ddatg) ; 

m_color(m_red) ; 

drawexp (dindatg, npts) ; 

end; 


{ 

procedure wrtfunc (wfunca : func_arr) ; 
var i : integer; 
begin 

writeln; writeCf ■ '); 
i 1; 

if wfunca [1] .datum * ends then writeln( 
while wfunca [i] .datum <> ends do 
begin 

with wfunca[i] do 
case datum of 

vars : write Cx’); 
coef : write (’ (C’ ,vali:i , ’) 
cons : write (* (A* ,vali: 1 , ') 
oper : case valo of 

oadd : write (’+’) 
osub : write(’-’) 
odiv : write(’/’) 
omul : write ( ’ * 1 ) 
opwr : write(’~’) 
otan : write ( ’tan’) ; 
osnh : write (’arctan’ 
end; 

end; 


- > 

{ write function in reverse } 

{ polish notation > 


undefined’) ;{ undefined functions > 

{ changing into human > 

{ readable forms > 


); 

); 


oexp 

write ( ’exp’) ; 

olog 

write ( ’log’) ; 

ocos 

write ( ’cos’) ; 

osin 

write ( ’sin’) ; 

oatn 

write ( ’arctan’ ) 

ocsh 

write ( 'cosh* ) ; 

otnh 

write('tanh’) ; 
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i :■ i + 1; 
end; 
end; 


{ > 

procedure readfunc(var rfunc : func_arr;var rcfterms , rconnum : integer); 
var ch : char; { read function in R.P.N. > 

i : integer; notaset : charset; 

begin 
page ; 

ref terms := 0; { reset function counts } 

rconnum := 0; 

i := 0; { initializations and } 

rfunc [1] . datum := ends; { instructions > 

notaset : * '-','*','/','" , ,*e','l', 

's' , 'c* , 'a* .’t'.’C’, ’S', 'T'.'x' 

writelnCto enter a constant, type fl then the constant # and <RETURN> * ) ; 
writelnCto enter a coefficient, type # then the coeff.# and <RETURN>*); 
writeln( ' other accepted notations 


writeln( * + , -, *, /, e(xp) , l(og), s(in) , c(os), x(var) , .(undo), ;(end‘); 
writeln( * t (an) , C(osh) , S(inh) , T(anh)'); 

writeln; writeln( ’Useful equations for statistical study’); 

writelnC Gaussian = (Cl) (C2) * (C3)/(2) (0 . 5) V (-0 . 5)x(C4) - (C3)/(2) "*exp* * ) ; 

writeln ( * Lorentz = (Cl) (C2) (2)/(2) ~*x(C3) - (2) *x(C2) (2)/- (2) A +/ ’ ) ; 

writeln( ’Resonance= (Cl) (2) *x(2) * (Cl) (2) (2) A x(C2) *(2) * + (0 . 5) V ’ ) ; 

writeln; writeln( 'enter function in REVERSE POLISH NOTATION’); 

writeln; write ('f - '); 

repeat 

read(ch) ; 

if not (ch in notaset) then { input error > 

begin 

write (' input error '); 
wrtfunc (rfunc) ; 

end 

else { interpreting inputs } 

begin 

i :* i + 1; 

rfunc [i+1] .datum := ends; 
with rfunc [i] do 
case ch of 

*x’ : datum := vars; 

: begin 

datum := cons; 
write ( ’ ? * ) ; 
readln(vali) ; 

if vali > rconnum then rconnum := vali ; 
wrtfunc (rfunc) ; 
end; 


esio 

February 14, 1987 


- 62 - 


’#’ : begin 

datum :* coef; 
write (’?’); 
readln(vali) ; 

if vali > rclterms then ref terms :* vali; 
wrtf unc (rfunc) ; 
end; 


i :* i - 2; 







begin 

datum : * 
case ch 

-\’e* 

oper; 

of 

•1* . ’s’ . ’c’ 

.’a’.’t’ 

9 C' ,'S’ , 

> ^ i 

’ + ’ : 

valo 

s 

oadd; 

’e* 

valo : = 

oexp; 

: 

valo 

r 

osnb ; 

’1* 

valo : = 

olog; 


valo 

s 

omul ; 

’c’ 

valo := 

ocos ; 

’/’ ! 

valo 

s 

odiv; 

’s' 

valo := 

osin; 

$ A | m 

valo 

m 

opwr ; 

•a* 

valo : s 

oatn; 

*t* : 

valo 

X 

otan; 

’C’ 

valo :* 

oesh; 

’S’ : 
end; 

valo 

s 

osnh; 


valo : = 

otnh; 


end; 

* ; ’ : datum : * ends ; 

otherwise ; 

end; 

end; 

until ch * ’ ; * ; 
end; 


{ > 

procedure current! it(x2 : real); { display infomation of fit } 


var i : integer; 

begin 

page; 

write (’Data file : ’ .name) ; 
case mode of 

stat : writelnC Sigma Y : Statistical’); 
none : writelnC No sigma y’); 
instr: writelnC Sigma Y : Instrumental’); 
end; 

write (’Function of graph is ’); { write function, its > 

wrtf unc (f unc a) ; { constants and coefficients } 

write In; 

for i := 1 to connum do 

writelnC A' ,i:l, * * *,cona[i]); 
for i :« 1 to cf terms do 

writeln(’C’ ,i:i, * ■ ’ .coeff [i] .a, ‘ Sigma C’ ,i: 1 , ’ = ’ .coeff [i] . siga) ; 
writelnC Chi square = ’ ,x2) ; 

writ e In ( ’Minimum y * ’.miny,* Maximum y * ’.maxy); 
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Maximum x = ’.maxx); 


writeln( 'Minimum x * ’.minx,' 
write(’hit any key to return '); 
read(mch) ; 
end; 

{ > 

procedure fit(var x2 : real); { controls fitting } 

var m : integer; xl.chidif : real; wt : oneddata; 

begin 

write(’color? ’); { input conditions for fitting} 

readln(color) ; 

write (’cut off point? ’); 

readlnCcutoff ) ; 

write(’enter maximal absolute difference of successive chi square (0.0001) ’); 
readln(chidif ) ; 

write ( ‘display graph as function being fitted? ’); 

readln(mch) ; 

lamda :* 0.001; 

nfree :* npts - cf terms; 

for m :■ 1 to cfterms do { initially guessed coeff. } 

begin 

write (’C* ,m:l, *?. delta? ’) ; 
readln(coeff [m] .a, coeff [m] . dela) ; 
end; 

for m:» 1 to npts do wt[m] :■ weight (indat p m) ; 

x2 :■ chisqr (indat .funca, coeff ) ; 

writeln(*chi square » ’ ,x2) ; 

if mch ■ ’y* then draweverything ( indat ,dat) ; 

repeat 

writelnC lamda ■ lamda,’ ....fitting....’); 

xl :■ x2; 

if nfree > 0 then curf it (indat .coeff ,x2 ,wt) ; 
writeln; writeln; writeln; 

for m := 1 to cfterms do { display temporary coeff. } 

writeln( ’ C ’ ,m: 1 , ’ * ’ , coeff [m] .a, 

’ Sigma C’,m:l,’ » ‘ .coeff [m] . siga) ; 
writeln(‘chi square m ’ ,x2) ; 
ifMnch * *y* then draweverything (indat ,dat) ; 
until abs(xl-x2) <* chidif ; •( stop fitting? } 

writeln( ’fitting completed’); 
draweverything ( indat ,dat) ; 


currentf it(x2) ; 
end; 

{ > 

procedure drawf unc ; { draw graph between xl & x2 } 

var m : integer; 
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begin 

m.clear; 

write (’enter minimum x and maximum x ’); 
readln(minx.maxx) ; 

dmaxy :* -9.999e99; { normalizing... > 

dminy :* 9.999e99; 

write(’New coefficients? ’); 

readln(mch) ; 

if mch = *y* then 

for m :* 1 to cf terms do 
begin 

write( ’C’ ,m: 1 , ' * *); 
readln(coeff [m] . a) ; 
end; 

corrx :■ long/(maxx - minx); 

dat :■ theodata(funca,coeff ,cona) ; 

maxmin; 

datg :* theonorm(dat) ; 

write( ’color of graph? ’); { drawing... > 

readln(color) ; 

drawf rame ; 

m.color (color) ; 

drawtheo(datg) ; 

end; 

{ - > 


procedure displaydata(dindat : expdata ; dnpts : integer) ; 


var i : integer; 

begin 

page; 

writeln( ’Data file: * .name) ; 
if dnpts > 0 then 

for i :■ i to dnpts do 
with dindat [i] do 
begin 

write(’X(* ,i:l, ’) - \x.’ 
write ( *Y( ’ ,i: 1 , *) - \y.’ 
writeln(’SIGMA-Y(’,i:l, ’) - 
if i mod 20 * 0 then 
begin 

write (’hit any key to 
read(mch) ; 
page; 
end; 

end; 

write(’hit any key to return’); 

read (mch) ; 

end; 


{ display current data file > 


’); 

’): 

’ ,sigy) ; 

{ end of screen } 


continue ’); 


{ end of data } 
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{ - - - > 

procedure drawdata(dindat :expdata) ; { draw data alone > 

var dindatg : expdata; 

begin 

maxmindata(dindat) ; { normalizing. . . > 

dindatg : = expnorm(dindat) ; 

m_clear; { drawing... } 

drawl rame ; 

m_color(m_red) ; 

drawexp (dindatg, npts) ; 

end; 

- — - > 


procedure dellunc(var dlunca : lunc_arr;var dclterms .dconnum: integer) ; 
var m : integer; 

begin { input a function > 

repeat 

readfunc (dlunca, dclterms, dconnum) ; { read function > 

wrtfunc (dlunca) ; 

writeln; 

write ( ’correct? ’ ) ; 
readln(mch) ; 

dlunca[0] .datum : = cons; 
dlunca [0] .vali :* 0; 

if (dconnum > 0) and (mch = *y # ) then { get constants il appropriate} 

lor m :* 1 to dconnum do 
begin 

write ( ’A* ,m: 1 , ** * ) ; 
readln(cona[m] ) ; 
end; 

until mch = 'y* ; 
end; 

{- - - } 


procedure transler(var tindat 

var i : integer; 

begin 

lor i :* 1 to npts do 
case dir ol 

1 : case subrec ol 

1 : tarray[i] :» 

2 : tarray[i] := 

3 : tarray[i] :» 
end; 

2 : case subrec ol 


expdata; var tarray: oneddata ; dir .subrec : integer) ; 

{ subprocedure of modify } 

{ procedure to transfer data } 
{ from one array to another } 

tindat [i] .x; 
tindat [i] .y ; 
tindat [i] .sigy; 
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1 : tindat[i].x :* tarray[i] ; 

2 : tindat [i] .y :* tarray[i]; 

3 : tindat [i] .sigy :■ tarray[i] ; 
end; 

end; 

end; 


procedure modify (var mindat : expdata) 
var i.tcom : integer; 

uarray , varray : oneddata; 
begin 
repeat 
page; 

writeln(’l --> 
writeln(’3 --> 
writeln(’5 --> 
writeln(’7 --> 
writeln(’0 --> 


(X -> Uarray) 

(Y -> Uarray) 

(SIGMA Y -> Uarray) 
(Uarray -> X) 

(Uarray -> Y) 

writeln( ’ ll--> (Uarray -> SIGMA Y) 
writeln( * 13--> (U - f(U)) 
writeln( ’ 15- -> (U * i(V)) 
writeln( ’other numbers--> end’); 
readln(tcom) ; 
case tcom of 

1 : transfer (mindat .uarray , 1 , 1) 

2 : transfer (mindat, varray. 1,1) 

3 : transfer (mindat .uarray , 1 ,2) 

4 : transfer (mindat .varray , 1 ,2) 

5 : transfer (mindat .uarray , 1 ,3) 

6 : transfer (mindat .varray , 1 ,3) 

7 : transfer (mindat .uarray, 2,1) 

8 : transfer (mindat, varray, 2,1) 

9 : transfer(mindat,uarray,2,2) 

10 : transfer(mindat.varray,2,2) 

11 : transfer (mindat .uarray, 2,3) 

12 : transfer (mindat .varray ,2,3) 
13,14,15: begin 

deff unc (funca ,cf terms .connum) 
for i :« 1 to npts do 
case tcom of 

13 : uarray [i] 


•( make modifications to data 
{ as a whole group 


{ instructions 


> 

> 


(X -> Varray) ’ ) ; 

(Y -> Varray) ’ ) ; 

(SIGMA Y -> Varray)’); 
(Varray -> X)’); 
(Varray -> Y) ’ ) ; 
(Varray -> SIGMA Y)’); 
(V = f (V)) ’) ; 


2 --> 

4 — > 

6 --> 

8 --> 

10 --> 

12 — > 

14— > 

16— > Display data’); 


{ transfering. . . 


14 : varray [i] 

15 : uarray [i] 
end; 

end; 

16 : displaydata (mindat ,npts) ; 
otherwise ; 


f (funca 
f (funca 
f (funca 


{ transforming data 


.uarray [i] , -le300 , coef f , cona) ; 
,varray[i] , -le300, coef f ,cona) ; 
, varray [i] ,-le300,coeff ,cona) ; 
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end; 

until (tcora < 1) or (tcom > 16); 
end; 

{ > 

procedure update(var uindat : expdata;var unpts : integer); 

var i,j ,k: integer; newval : real; { make changes and add data > 

begin 

repeat 

writeln( ' 1-Remove data point 2-Add data 3-Change data 4-Display data ' ); 
write( 'other numbers to exit ’); { updating options > 

readln(i) ; 
case i of 

1 : repeat { removing a single data point} 

write ( 'remove item #? (9999 to terminate) ' ); 
readln(j) ; 

if (j<unpts) and (j>0) then 
begin 

for k := j to unpts -1 do 
uindat [k] : = uindat [k+1] ; 
unpts unpts - 1; 

end 

else if j - unpts then unpts : = unpts - 1; 
until j * termx; 

2 : readdata (uindat .unpts , unpts) ; { add data points } 

3 : repeat { change data } 

write ('item #? (9999 to terminate) ’); 
readln( j ) ; 
if j <= unpts then 
with uindat [j] do 
begin 

write( 'X( * , j : 1 , ' ) = '.x.' '); 

write ( 'Y( * , j : 1 , *) = \y.’ *); 

writeln(*SIGMA-Y(' , j :1, *) = \sigy); 
write( 'Change 1-X 2-Y 3-SIGMA-Y ? *); 

readln(k) ; 

write ('new value? '); 
readln(newval) ; 
case k of 

1 : x newval; 

2 : y : = newval ; 

3 : sigy := newval; 

otherwise; 

end; 

end; 

until j * termx; 

4 : displaydata (uindat .unpts) ; { display data } 

otherwise ; 
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end; 

until (i<l) or (i>4) ; 
end; 

i -> 

function diffdata(dindat : expdata ; ddat : theograf ) :expdata; 

var i : integer; { get difference plot > 

begin 

diffdata :* dindat; 
for i :* 1 to npts do 

diffdatati] .y :■ dindat [i] .y - ddat [round( (dindat [i] .x-minx)*corrx + left)]; 
end; 


{ > 

procedure getfunc(var gfunc : func_arr;var gcf terms, gconnum : integer); 
var i : integer; { read function in R.P.N. } 

funcname : string [12] ; dumfunc : filefuncrec; { from file > 


begin 

page; 

write('Read function from file? ’); 
readln (funcname) ; 
reset (fileb, funcname) ; 

gcf terms :* 0; { initializations > 

gconnum :* 0; 
i :* 0; 

gfunc [0] . datum :* cons; 
gfunc [0] . vali :* 0; 
repeat 

i :* i + 1 ; 

gfunc [i+i] .datum :■ ends; 
dumfunc :■ fileb*; 
with gfunc [i] do 

with dumfunc do { intepreting saved function > 

case sort of 

1 : datum :* vars; 

2 : begin 

datum :* cons; 
vali :■ sub; 

if vali > gconnum then gconnum : = vali; 
end; 

3 : begin 

datum :* coef ; 
vali : * sub ; 

if vali > gcf terms then gcf terms :* vali; 
end; 

4 : begin 

datum :■ oper; 


- 69 - 


CS10 

Ffebruary 14, 1987 


case sub of 

1 : valo := oadd; 

2 : valo := osub; 

3 : valo :» omul; 

4 : valo := odiv; 

5 ; valo := opwr; 
10: valo := otan; 
13: valo := osnh, 

end; 

end; 

otherwise ; 
end; 

get(fileb) ; 
until eof (f ileb) ; 
close(f ileb) ; 
wrtf unc (gf unc ) ; writeln ; 
if gconnum > 0 then 

for i := 1 to gconnum do 
begin 

write ( ’ A ’ , i : 1 , ' = '); 
readln(cona[i] ) ; 
end; 

end; 

{ 


6 : 

valo 

: = 

oexp; 

7 : 

valo 

: = 

olog; 

8 : 

valo 

: = 

ocos; 

9 : 

valo 

: = 

osin; 

11: 

valo 

: = 

oatn; 

12: 

valo 

: = 

ocsh; 

14: 

valo 

: = 

otnh; 


{ get constants if appropriate} 


procedure savefunc (sfunc : func_arr) ; { save function in R.P.N. } 

var i : integer; 

funcname : string [12]; dumfunc : filefuncrec; 

begin 
page; 

write ('Save function to file? '); 
readln(funcname) ; 
rewrite (f ileb .funcname) ; 


i := 1; 
repeat 

with sfunc [i] do 
case datum of 


vars : 

dumfunc . sort : = 

l; 


cons : 

begin 




dumfunc . sort 

: = 

2; 


dumfunc . sub 

* s 

vali ; 


end; 



coef : 

begin 




dumfunc . sort 

: = 

3; 


dumfunc . sub 

; s 

vali ; 


end; 



oper : 

begin 




dumfunc . sort 

; s 

4; 


{change function into textfile} 
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with dumfunc do 
case valo of 


oadd 

sub 

; s 

i 


oexp 

sub 

: * 

6; 

osub 

sub 

:« 

2 


olog 

sub 

* s 

7; 

omul 

sub 

: = 

3 


ocos 

sub 

: = 

8; 

odiv 

sub 

: = 

4 


osin 

sub 

* x 

9; 

opwr 

sub 

| X 

5 


oatn 

sub 

• X 

11; 

otan 

sub 

♦ X 

10; 

ocsh 

sub 

: = 

12; 

osnh 

sub 

:* 

13; 

otnh 

sub 

:» 

14; 


end; 

end; 

otherwise ; 
end; 

fileb* : = dumfunc; 
put(f ileb) ; 
i :« i + 1; 

until sfunc [i] . datum * ends; 

close (f ileb, ’save’) ; 

end; 

{ — - - > 

procedure extrapolate; { extrapolate known function > 

var ex,ey : real; 

begin 

page; 

writeln(’To terminate, enter 9999 for x’); 
repeat 

write ( 'x * ’ ) ; 
readln(ex) ; 

ey : = f (funca, ex, cutoff .coeff ,cona) ; 
writeln(’x = *,ex,’ f(x) = ’,ey); 
until ex = termx; 
end; 

{ - - — - > 

procedure wait_for_pen(var wpen : m_tablet_info; wcorx.wcory ,wx,wy : real); 
begin { wait til pen is depressed > 

repeat { and return coordinate > 

repeat 

m_readpen(wpen) ; 
until wpen. moving or wpen.dn; 
page; 

writeln(wpen.x*wcorx + wx, ’ ’ ,wpen.y*wcory + wy) ; 

until wpen.dn; 
end; 
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procedure readgraf (var rindat : expdata;var rnpts : integer); 

var mypen : m_tablet_inf o ; { read data from graph > 

rmaxx, rmaxy, rminx, rminy, rcorx, rcory : real; { on graphics tablet > 
amaxx, amaxy, aminx, aminy, scorx, scory : real; 

begin 

page; 

m_move(10 ,360) ; 

m_displaytext( 'place graph on graphics tablet*); 
m_move(10 ,350) ; 

m_di splay text ( *move pen to the lower left corner of your graph and press*); 

wait_f or_pen (mypen, 1 ,1,0,0); 

rminx := mypen. x; 

rminy := mypen. y; 

m^clear; 

m_move(10 ,360) ; 

m.displaytext ( 'move pen to the upper right corner of your graph and press*); 
wait.f or _pen (mypen, 1 ,1,0,0); 
rmaxx :* mypen. x; 
rmaxy := mypen. y; 

write(*enter actual minx, maxx, miny, maxy *); { get actual values on graph > 

readln(aminx, amaxx .aminy, amaxy) ; 

rcorx := (amaxx - aminx)/ (rmaxx - rminx); { normalizing... } 

rcory := (amaxy - aminy) /(rmaxy - rminy); 

scorx :* long/(rmaxx - rminx); 

scory := height/(rmaxy - rminy); 

m_clear ; 

mjmove(10 ,360) ; 

m_displaytext ( 'move to point to be collected and depress, box 1 to end*); 

rnpts : - 0; 

m_color(m_yellow) ; 

m_drawrect(lef t .bottom, right .top) ; 

m_color(m_red) ; 

mch :* *n*; { flag to stop data entry > 

repeat { collecting data points > 

rnpts :« rnpts + 1; 
wait _f or_pen (mypen , rc orx , rc ory , 

-rminx*rcorx + aminx, -rminy*rc ory + aminy); 
rindat [rnpts] .x :■ ( mypen. x-rminx)*rcorx + aminx; 

rindat [rnpts] .y :• (mypen. y-rminy)*rc ory + aminy; 

rindat [rnpts] . sigy :■ 0; 

m_circle (round(scorx* (mypen. x-rminx)+left) , 

round(scory* (mypen. y-rminy)+bottom) ,1) ; 
if ( (mypen. x>0) and (mypen. x<22) ) and 

( (mypen. y>393) and (mypen. y<417) ) then begin 
write('data entry complete? *); 
readln(mch) ; 
end; 

until (mch 88 *y*) or (rnpts * 100); 
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rnpts :■ rnpts - 1; {do not enter menu value > 

end; 

{ - - > 


2-Save current file 
5-Fit 

8-Display data 
11-Difference Plot 
14-Save function 
17-Digitize graph 


proc edure menu ; { 

var com : integer; difdat : expdata; { 

begin 

writeln( ’ 1-Get a file 
writeln( '4-Input data 
writeln( '7-Plot function 
writeln( ’ 10-Update data 
writeln( ’ 13-Get function 
writeln( ’ 16-Chi square 
writeln( ’ 19-Exit program’); 
readln(com) ; 

case com of { 

1 : begin { 

readf ile(indat f npts,0) ; 
nfree :» npts - cf terms; 
end; 

2 ; if npts > 0 then savef ile(indat ,npts) ; { 

3 : deffunc (funca.cf terms, connum) ; { 

4 : begin { 

name : * * not named ’ ; 
readdata(indat ,npts ,0) ; 
nfree : • npts - cf terms; 
end; 

5 : if npts > 0 then fit(x2); { 

6 ; if npts > 1 then drawdata(indat) ; { 

7 : drawfunc ; { 

8 : displaydata(indat ,npts) ; { 

9 : if npts > 0 then modif y(indat) ; { 

10: if npts > 0 then update(indat.npts) ; { 

11: begin { 

writeln( 'enter c for difference plot IFF 

readln(mch) ; 
if mch ■ 'c' then 
begin 


difdat :» diffdata(indat.dat) ; 
maxmindat a (difdat) ; 
if npts > 1 then drawdata (difdat) ; 
currentf it(x2) ; 
end; 


display menu and 
interact with user 

3-Define a function’); 
6-Plot data’); 

9-Modify data’); 
12-Function k Data’); 
15-Extra- ( Inter) -polate ’ ) ; 
18-Combine files’); 


executing chosen option 
get data file 


save datafile 
user input function 
user input data 


fitting 
plotting data 
plotting function 
display data 
modify group of data 
update one data point 
difference plot 
data was fitted’); 


end; 

12: begin { plot both data and graph 

writeln( ’enter c for plot IFF data was fitted’); 
readln(mch) ; 

x2 :* chisqr(indat,funca,coeff) ; 


> 

> 


> 

> 


> 

> 

> 


> 

> 

> 

> 

> 

> 

> 


> 
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if inch * 'c* then begin draweverything ( indat , dat) ; 
current f it (x2) ; end; 

end; 

13: getfunc (funca.cf terms , connum) ; { read function from file > 

14: if funca[l] . datum <> ends then savefunc (funca) ; { save defined function > 
15: extrapolate; { extrapolating from known fen} 

16: begin { given a function and a set } 

draweverything(indat ,dat) ; { of data point, plot and } 

x2 := chisqr(indat .funca.coeff ) ; { determine chi square } 

currentf it (x2) ; 
end; 

17: readgraf (indat ,npts) ; { read from graph on tablet } 

18: if npts > 0 then readfile (indat ,npts,npts) ;{ Combine two files } 

19: exitprog := true; { exiting program } 

otherwise ; 
end; 
end; 

{ } 


procedure currentinfo; { display status of data etc. } 

var i : integer; 
begin 
page ; 

write ( ’Current data file is * .name) ; { name of data and function } 

case mode of 

stat : writeln(* Sigma Y : Statistical'); 
none : writelnO No sigma y'); 
instr: writelnC Sigma Y : Instrumental'); 
end; 

write ( 'Current function is ’); 

wrtfunc (funca) ; 

writeln; 


if connum > 0 then for i := 1 to connum do { constants } 

writeln('A* ,i:l, ’ = ’,cona[i]); 
writeln; 
end; 

{ > 

begin 

page; { initialization } 


funcatl] .datum := ends; 
name := 'none'; 
exitprog :* false; 
cf terms := 2; 
mode : = none ; 
m_init_graphics ; 

repeat { begin interaction with user } 
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current info; 
menu; 

until exit prog; 
end. 
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